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A new type of asymptotic behavior in a game dynamics system is discovered. The system exhibits 
behavior which combines chaotic motion and attraction to heteroclinic cycles; the trajectory visits 
several unstable stationary states repeatedly with an irregular order, and the typical length of the 
stay near the steady states grows exponentially with the number of visits. The dynamics underlying 
this irregular motion is analyzed by introducing a dynamically rescaled time variable, and its relation 
to the low- dimensional chaotic dynamics is thus uncovered. The relation of this irregular motion 

■ with a strange type of instability of heteroclinic cycles is also examined. 
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In this paper, we will report a new type of asymptotic behavior observed in a game dynamics system. The behavior 
exhibits a combined nature of the attraction to a heteroclinic cycle and a chaotic motion; the trajectory visits some 
[T i ■ of the saddle points with an irregular order and the lengths of the stay near the saddles are also irregular. The 
\ irregularity here is related to an intrinsic chaotic dynamics and not due to the randomness from noise. The behavior 

■ is robust against a small variation of the parameters, and has a close relation with the structurally stable heteroclinic 
cycles studied recently in various contexts jlj Q || [Q . 

The game dynamics system we will study is a kind of population dynamics, and is represented with a set of equations 

> 

in; d_ 

S : *i = Sl - s, (2) 

o ■ 

s i = 9ij x ji (3) 

ON ■ i=i 



s = J2 s * x *> ( 4 ) 

-a 
i 

O ' with a constraint. 

O- 5^3* = 1. (0<Xi<l). (5) 
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Xi represents the renormalized population of species i, s, is the score (or fitness) of species i, and s is the average 
score over the population, {gij} defines the interaction between species, each representing the score of species i in 
. . . i its battle with j. 

Such a model was proposed by Taylor and Jonker in the context of the evolution of strategy into ESS (evolutionary 
stable states) || ; it is also regarded as a simplified model equation for the molecular evolution ^ , and gives a minimal 
model for the system of interacting self replicators . It is also related to the Lotka-Volterra equation through a 
transformation of variables |t]] . 

The system exhibits various types of asymptotic behavior. The simplest one is the relaxation to a resting state. 
The chaotic oscillation as well as regular oscillation can be observed in some range of parameters when 4 or more 
species are involved || ^ [jlO| . Attractive heteroclinic cycles also are observed generically for systems with 3 or more 
species [fill . 

As is easily seen from the constraint (|^) , the phase space of the game dynamics system with n species is given by 
an (n — l)-dimensional simplex. The surface of the simplex consists of n hyper-planes Xi = 0, each of which gives an 
invariant set. Therefore, any intersection of these hyper-planes is also an invariant set. The flow on such an invariant 
set is described by the game dynamics equations of the original form but with a reduced number of species, and gives 
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a subsystem of the original one. Note that there exist many ( at least n ) fixed points on the surface of the simplex, 
since each of these subsystems has at least one fixed point. 

A peculiar feature of this system is the existence of robust heteroclinic orbits. This is basically due to the hierarchical 
structure of the invariant set of the system described above. Suppose that two saddles A,B and a heteroclinic orbit 
from A to B exist on an invariant set which corresponds to a subsystem with d species. The heteroclinic orbit will be 
robust against a small variation of the parameters if 

dim A" +dimB s > (d- 1), (6) 

where A u denotes the unstable manifold of saddle A in the subsystem [[12|| . Thus, if the saddle B is an attractor in 
the subsystem, i.e., dim£> s = [d — 1), the heteroclinic orbit is robust. We will concentrate on such robust heteroclinic 
orbits, therefore, every saddle which appears in this paper is an attractor in the subsystem consisting of the species 
with non-zero Xi at the saddle point. Note that such saddles and robust heteroclinic orbits are all on the border of 
the phase space, because it must belong to a certain subsystem. 

Some of the results presented in this paper are obtained from numerical simulations. A small noise in the calculation 
may drastically change the behavior of the trajectory especially in the vicinity of the heteroclinic orbits. To suppress 
the numerical error, the simulations were done with algorithm which keeps a constant precision for them. It is realized 
by doing the numerical simulations of the dynamics of loga^ which is, in principle, equivalent to the original dynamics 
of Xi. The error for the population near the zero value is drastically reduced by keeping tracks its logarithm. In 
the contrast, the error for non-zero valued population is automatically reduced by the dynamical contraction in the 
relevant direction, as far as we work with robust heteroclinic orbits mentioned above. 

We will now look into the system behavior. An example is presented in fig.l. It is obtained from the simulation for 
a system with 5 species with interaction matrix 
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where X = 0.32. The population are set as (1/5, 1/5, 1/5, 1/5, 1/5) initially. The behavior looks like the attraction 
to an heteroclinic cycle. The trajectory visits several quasi-stable states, and the period of the stay there shows a 
geometrical expansion. There is, however, irregularity in the order of visits and lengths of the stays, which is a major 
difference from typical cases of the attraction to an heteroclinic cycle (fig.||) . The transition between these two types 
of behavior, namely the irregular 'rambling' motion and the regular attraction to heteroclinic cycle, will be studied 
afterwards using the same matrix as the above but with different value of X . 

We will analyze the asymptotic dynamics of this system. Here we assume that the trajectory stays sufficiently 
long in the neighborhoods of the saddles, while the transitions between them take place in a much shorter time. The 
irregular rambling motion and also the attraction to a heteroclinic cycle are analyzed as follows. In the analysis, we 
will work with the logarithm of the populations 

yi= log Xi, (8) 
and consider an approximation for the dynamics of yi. Here we write the approximated dynamics in the form 

jY, = MY), (9) 

where Yi is used instead of yi to avoid confusion and Y denotes (Yi, Y~2, ••• ) ^n)- It is assumed to be valid in the 
vicinity of saddles and/or robust heteroclinic cycles. All of those relevant saddles and heteroclinic orbits are on the 
border of the phase space by assumption. Therefore, some of the species have nearly zero valued population in the 
considering region. We will call such species as 'inactive' species and the others as 'active' ones. 
The considering simplification consists of following two points, namely, 

• ignoring the change of Yi made during the transition between saddles, 

• linearizing the flow in the neighborhood of saddles. 

The dominant contribution to the variation of yi comes from the stays in the neighborhoods of saddles, if the transi- 
tion between saddles is sufficiently swift as assumed. Thus we will take an approximation that the transition between 
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saddles is an instantaneous event, and neglect the change of Yi during the transition. With such a simplification, 
Yi contains non vanishing error inevitably, and it is useless to consider a direct correspondence between Yi and Xi. 
Thus we simply assume that Y{ is zero when species i is active, and Yi < corresponds to a;, » 0. In this case, the 
constraint on Y corresponding to (|5|) is represented as, 

maxY- = 0. (10) 

i 

As is mentioned above, the transition between saddles is neglected in the dynamics of Y. Therefore, all of the change 
of Y is made during the stay in the neighborhoods of certain saddles. By linearizing the flow in the neighborhood of 
saddles, / can be written formally as, 

MY) = \f, (11) 

where A denotes the neighboring saddle which should be determined from Y, and \f is the growth rate of species i 
at saddle A. 

Here we define / by simply assuming that A is determined as the attractor of the subsystem which consists of 
the species with Yi — 0. Thus invasion of species I corresponds to an instantaneous event of Yj reaching 0, and 
immediately after that the system is assumed to be in the neighborhood of a new saddle. Note that this definition 
of / may be invalid for some values of Y, where two or more attractor coexist in the associated subsystem. This 
definition, however, is still valid for most value of Y and it is sufficient for the following argument. 

Although \f is zero if xf has nonzero value, Yi — does not necessarily lead to Y — 0. This is because the 
attractor can be on the border of the phase space of the associated subsystem, that is, xf can be zero even if Yi is 
zero. In any case, the constraint ( |Tl| ) is not violated because Xf cannot be positive if Yi is zero, in order for the fixed 
point A to be an attractor in the associated subsystem. 

Since the function f(Y) depends only on whether each component of Y is zero or not, it can be formally expressed 
as a function of ff = Y/\Y\. Thus, by using the notations 

ns w (12) 

U=\Y\, (13) 

m = M-{ff'M)n, (14) 

h(ff) = ff-f(ff), (15) 

a set of equations, 

j t fi=u- i m, (is) 

±U = h(f f) , (17) 



is obtained. It should be noted that, corresponding to fll 1|) , there is a constraint on ff, such that, 

m&xrii = 0. (18) 

i 

We now introduce a dynamically rescaled time variable r defined as 



Then the equations (|16D,(|17j) can be rewritten as 

£n=M, (20) 

^-(logU) = h(ff). (21) 

Thus we obtain an autonomous equation for ff, and this determines the sequence of the saddles visited in the rescaled 
time t. On the other hand, the time scale of the motion which associated with the length of the stay near saddles is 
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determined by U. As can be seen from ( |20| ) and (pl|), the motion has no characteristic time scale. Steady motion of fj 
in r results in a drift of (log U) with some average velocity measured with r. This results in a geometrical change of 
the time scale of the motion in t through ([l9|) . Thus a limit cycle solution of (|2(]) corresponds to a heteroclinic cycle 
in the original game dynamics system. 

So far, it has not been questioned whether the trajectory of rf(r) is regular or chaotic. We now look into this point 
more carefully with some case studies. For this purpose, it is convenient to work with discrete maps which represent 
the change of Y in a neighborhood of a saddle, rather than the continuous flow of Y or if. 

The change of Y near a saddle is represented by a piecewise-linear map. This is because when the invader is 
specified as 7-th species, the map is expressed in a linear form, 

y(") , > y("+!) _ y{ n ) _ ^i_y( n ) (22) 

Xf 

where A : I denotes the saddle and invading species respectively, and fW denotes the value of Y at the n-th transition. 
The invader is determined in accordance with the value of Y^ n \ thus the map becomes 'piecewise'-linear as a whole. 

The return map for Y at some Poincare section can also be represented by a piecewise-linear map. Actually, if the 
sequence of saddles — or equivalently the sequence of invaders — is specified, the return map for Y is obtained as a 
product of the linear maps with the same form as (p2[). Thus each segment of the return map is a linear map. The 
sequence of the saddle is determined in accordance with the value of Y at the cross section, and the return map also 
becomes piecewise-linear. 

Let us consider the specific system (Q) and look into its asymptotic behavior. A Poincare section at the transition 
path from saddle {3,4} (where species 3 and 4 is steadily coexisting) to saddle {5} is considered. The saddle {3,4} 
is apart from the hyperplanes X3 — and 2; 4 = 0, and after the invasion of species 5, the orbit gets apart also from 
X5 = 0. Therefore, the value of Y at the corresponding moment has only two non-zero components, Y\ and Y%- As 
is obvious from the above argument, the sequence of the saddles visited remains unchanged if Y is multiplied by a 
scalar factor. Thus the sequence of the saddles can be determined from the ratio between these components, i.e., 

Z = Y 1 /Y 2 . (23) 

In this way, we can analyze the saddle-sequence with the return map for Z . If the chaotic trajectory is produced by 
this map, the corresponding saddle-sequence will also be chaotic. 

The return map for Y and also the return map for Z can be calculated analytically from the interaction matrix 
given in (Q). The map for Z is presented in fig.[| The part of the map displayed consists of three elements, each of 
which corresponds to a particular sequence of the saddles. The three types of the sequence are presented in table Q. 
The map for Z has one fixed point associated with heteroclinic cycle C2 which is an unstable fixed point, and the 
map generates a chaotic sequence of CI, C2 and C3 under general initial conditions. 

The sequence of the saddles is thus expected to become chaotic in this case, which is confirmed in the numerical 
simulation. A return map for z = 2/1/2/2 is studied, and its Poincare section is taken at X5 = e, £5 > 0, X3 > e, £4 > e 
with e = 10 -8 . This hyper-plane corresponds to the instant of the 'invasion' of species 5 to the saddle {3,4}. The 
behavior of this map is well reproduced by the return map for Z which can be calculated analytically, as is exhibited 
in fig.||. It is thus clear that the irregularity of the sequence of the saddles observed in figsjljand ^|is not due to noise 
or numerical errors but generated by the intrinsic chaotic dynamics. 

We have reduced the original 4 dimensional flow to a one dimensional map. One of the variables is eliminated by 
taking a Poincare section of the flow. The reduction of the second variable is enabled by the strong contraction of the 
flow, and we could separate the last one because of the linear approximation. Note that it depends on the particular 
system or particular cycle, that how many variables can be reduced by the contraction of the flow. If all the relevant 
heteroclinic orbits belong respectively to 1 dimensional subsystems, no variables can be reduced by the contraction. 

Finally we will analyze the break down of a heteroclinic cycle attractor, which leads to the emergence of the chaotic 
rambling motion. 

The return map for Z (fig.||) has a fixed point associated with heteroclinic cycle C2 as stated before. The stability of 
this fixed point is related to the property of the matrix M^ C2 ^ which represents the return map for (Yi, YV) associated 
with the same cycle. The fixed point of the return map for Z corresponds to an eigenvector of M^ C2 \ and the 
condition for the cycle C2 to be an attractor is expressed in terms of the eigenvalues and eigenvectors of M^ C2 \ that 
is, 

• A F > 1, 

. (VMF),|A F | > |A,|, 
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where Ap is the eigenvalue associated with the fixed point of the return map for Z. 

Even if the first condition is not satisfied, the fixed point in the return map remains stable. In this case, however, the 
iteration of this cycle makes \Y\ smaller and smaller, which means that the trajectory diverges from the heteroclinic 
cycle gradually. On the other hand, when the second condition is violated, the fixed point becomes unstable. Since 
the map for Y is linear in the region associated with C2, there is no nonlinear effect to suppress the instability of the 
fixed point, so that the trajectory will inevitably be thrown out of C2. In this way, the destabilization of the fixed 
point in the return map leads directly to the loss of attractiveness of heteroclinic cycle C2. 

Matrix M^ C2 ^ can be calculated from the matrix in (0) as, 

M (C2) _ ( 1-9389 + 6.418A 1.2386 - 12.836A \ , > 

~ ^ 2.841 + 8.42A 0.634 - 16.84A ) ' [ ' 

The fixed point of the return map for Z corresponds to the largest eigenvalue of this matrix which is larger than 1 for 
the considering range of X, so that cycle C2 is attractive if the fixed point of the map is stable. The absolute value 
of the negative eigenvalue, however, can exceed the largest eigenvalue, and destabilization of the fixed point can be 
observed by changing parameter X. Its critical value X c is obtained from the condition, 

Trace M (C2) = 0, (25) 

and we obtain X c = 0.24687 ■ ■ •. If X is larger than X c , the fixed point is unstable and the observed saddle-sequence 
becomes chaotic. On the contrary, the attraction to the heteroclinic cycle C2 is generic behavior when X < X c . The 
second return map for Z (fig.||) clearly shows that when the fixed point of the return map becomes unstable, the cycle 
C2 ceases to be an attractor and then a chaotic rambling motion appears. 

The stability of the heteroclinic cycle has rather peculiar features, as also noted by Melbourne Even in the 
case of X < X c , the heteroclinic cycle C2 is not stable in the usual sense, because an arbitrarily small additive noise 
can kick the trajectory out of C2 into a 'transverse' direction. It should also be noted that even in the case of X > X c , 
there exists a set of initial conditions with zero measure corresponding to the unstable fixed point in the map for Z, 
which is attracted to the cycle C2 (that is, the 'stable manifold' of C2). Another interesting phenomenon related to 
stability of the structurally stable heteroclinic cycles is reported recently . That is the coexistence of two or more 
heteroclinic cycle attractors sharing a common heteroclinic orbit. It will correspond to the case when the return map 
for ff has two or more stable fixed points. 

We have worked with the case of 5 species, because it is a minimal system capable of displaying an irregular saddle- 
sequence as an asymptotic behavior. Since ff has two constraints, \ff\ = 1 and max 7^ = 0, the dynamics of ff can be 

i 

represented by a flow on a n — 2 dimensional spherical surface (with n denoting the number of species). Therefore 
the asymptotic behavior of ff cannot be irregular if n is smaller than 5. It should be also noted that 5 species are 
required also in Lotka-Volterra systems for the emergence of the irregular rambling motions, if we assume that the 
total population always has finite value. 

We will mention the nature of the approximation used in this paper. The trajectory of ff gives a 'good' approximation 
for the trajectory of y, in the sense that for any positive integer k and for generic value of the initial condition ffo, 
there exist a set of initial conditions of y which give the saddle-sequences which have opening k saddles in common 
with the one associated with the trajectory of the ff. It can be proven with assuming the non-degeneracy of {gij}. A 
rough sketch of the proof is given below. The complete proof will be given in another paper fl5| . 

We can confirm the existence of an upper limit for error associated with the map ( p2[ ) (where the error means 
the difference between the change of Y made by this map and the change of y made during the corresponding stay 
and transition) which does not depend on the value of V. It is also confirmed that, as is obvious from (|22] ) and the 
assumption of non-degeneracy, an upper limit exists for the magnification rate of the deviation of two trajectories of Y 
in each step. Therefore, we can evaluate the upper limit of the difference between the value of y and its approximation 
Y after k transitions. With considering the fact that the map (]2^) is linear, we can choose the initial condition Yq, 
s.t. i) the corresponding saddle-sequence is the same as the one given by the trajectory started from 770, and ii) the 
opening k saddles of the associated saddle-sequence is not changed even if the error is taken into account. Y = U^ff 
with a sufficiently large Uq gives an example of the 'good' initial condition for Y, and the trajectory with the initial 
condition for y corresponding to the Yq gives the saddle-sequence with the desired property. 

As is mentioned above, we have analyzed the global behavior of the trajectory basically in the logarithmically 
transformed coordinate. Such a treatment is enabled by the robust existence of the hierarchically structured invariant 
sets (Yl^i — 0, which form the border of the phase space), and does not rely on a specific form of interaction. Thus 
the same method can be applied to the analysis of the system with nonlinear interaction as well. As is reported by 
several authors, structurally stable heteroclinic cycles can exist in systems with a certain kind of symmetry. Such 
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systems have invariant sets with similar structure to that in the game dynamics system. Thus the chaotic rambling 
over quasi stable states could also be observed in such systems. 

The relation between Y and / is not simple in the present system, and it is difficult to write down explicitly the 
rescaled dynamics (p0[). The dynamics of the transitions among saddles can be treated more elegantly in the case of 
a coupled logistic equations, where the dynamics with a rescaled time is written explicitly, and the emergence of a 
chaotic saddle-sequence is demonstrated Q . 
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TABLE I. Relation between the value of Z and the sequence of saddles. 
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FIG. 1. An example of the temporal evolution of the population of species, in the case of attractive heteroclinic network. 
The population of all 5 species are plotted, {gij} is chosen as in (Q) with X = 0.32, and the initial condition is set as 
xi — X2 — ■ ■ ■ = 0.2. Numerical simulation is done with the Runge-Kutta method with variable step size, taking loga;; as 
variables. 

FIG. 2. A temporal sequence of the saddles for a regular (a: A = 0.24), and an irregular motion (h:X — 0.32). The abscissa 
represents time in the log scale, and the ordinate displays the nearest saddle. The transition between saddles is omitted in this 
figure, however, the process is so relatively fast that it constitutes practically zero width with this scale. The initial condition 
is the same as in fig.l. 

FIG. 3. First return map for Z (solid line; analytically calculated), and z = 2/2/2/1 (plotted with circles; obtained from a 
numerical simulation). The parameter X is set as X = 0.32, and the considered Poincare section corresponds to the exit from 
saddle {3,4} towards saddle {5}. The section is taken as X5 = e, is > 0, (xs > e, X4 > e) with e = 10~ 8 in the numerical 
simulation. Plotted data corresponds to successive 200 turns (corresponding range of logi is [27.5, 105.25]). 

FIG. 4. The second return map for Z, obtained for A = 0.24 (dashed line) and 0.25 (solid line). The considered section is 
the same as in fig. 3. 
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